#directory, packages and functions
#set working directory
setwd("~/Documents/RSfPETCog")
#list packages
pags<-list("openxlsx", "reshape","tidyverse","ggpubr","rstatix", "plyr", "corrplot", "ggplot2","reshape2","lme4","gghalves")
# packages for: loading xlsx files
#long/wide tranform
#load packages
lapply(pags, require, character.only=TRUE)
## Loading required package: openxlsx
## Loading required package: reshape
## Loading required package: tidyverse
## Warning: package 'ggplot2' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks reshape::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::rename() masks reshape::rename()
## ✖ lubridate::stamp() masks reshape::stamp()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: ggpubr
##
## Loading required package: rstatix
##
##
## Attaching package: 'rstatix'
##
##
## The following object is masked from 'package:stats':
##
## filter
##
##
## Loading required package: plyr
##
## ------------------------------------------------------------------------------
##
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
##
## ------------------------------------------------------------------------------
##
##
## Attaching package: 'plyr'
##
##
## The following objects are masked from 'package:rstatix':
##
## desc, mutate
##
##
## The following object is masked from 'package:ggpubr':
##
## mutate
##
##
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
##
##
## The following object is masked from 'package:purrr':
##
## compact
##
##
## The following objects are masked from 'package:reshape':
##
## rename, round_any
##
##
## Loading required package: corrplot
##
## corrplot 0.92 loaded
##
## Loading required package: reshape2
##
##
## Attaching package: 'reshape2'
##
##
## The following object is masked from 'package:tidyr':
##
## smiths
##
##
## The following objects are masked from 'package:reshape':
##
## colsplit, melt, recast
##
##
## Loading required package: lme4
## Warning: package 'lme4' was built under R version 4.3.3
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## The following object is masked from 'package:reshape':
##
## expand
## Warning in check_dep_version(): ABI version mismatch:
## lme4 was built with Matrix ABI version 1
## Current Matrix ABI version is 0
## Please re-install lme4 from source or restore original 'Matrix' package
## Loading required package: gghalves
## [[1]]
## [1] TRUE
##
## [[2]]
## [1] TRUE
##
## [[3]]
## [1] TRUE
##
## [[4]]
## [1] TRUE
##
## [[5]]
## [1] TRUE
##
## [[6]]
## [1] TRUE
##
## [[7]]
## [1] TRUE
##
## [[8]]
## [1] TRUE
##
## [[9]]
## [1] TRUE
##
## [[10]]
## [1] TRUE
##
## [[11]]
## [1] TRUE
##Define functions
#function for returning individual cor.test pvalue per cor in data matrix compatible with corrplot
rcor.pvalue<-function (mat, ...) {
mat <- data.matrix(mat)
p <- ncol(mat)
index <- t(combn(p, 2))
nindex <- nrow(index)
p_mat<- matrix(nrow = p, ncol = p)
for (i in 1:nindex) {
p_mat[index[i,1], index[i,2]]<-cor.test(mat[, index[i,1]], mat[, index[i,2]], ...)$p.value
}
pMat <- p_mat
colnames(pMat)<-colnames(mat)
rownames(pMat)<-colnames(mat)
diag(pMat)<-0
pMat
}
#define stats parameters to be calculated by aggregate
myStats<-function( x){
c(Mean = mean ( x),
SD = sd( x))
}
#read data
CogNET<-read.xlsx(xlsxFile = "RSNNetworkROIsextracted.xlsx", sheet = "Tabelle1",
colNames = FALSE)
colnames(CogNET)<-c("MPFC","LPL","LPR","PCC","SML","SMR","SMS","ACC","AIL","AIR","PFCL","PFCR","SMGL","SMGR","LPFCL","PPCL","LPFCR","PPCR","CA","CP")
#read global mean values from file
my_globals <- read.delim("globals_3mm.txt",header=FALSE, sep=",")
#Intensity normalization, divide uptake in each ROI by global per timeframe
CogNETp<-CogNET/my_globals[1:810,]
#export seed values from MPFC after intensity-norm
#export<-CogNETp[,c(1:20)]
#write.table(export$MPFC, "Cog_Net_gm.txt", col.names = FALSE,row.names=FALSE)
#define subject variable subject id
subject_id<-c(rep("RS_fPET001",30), rep("RS_fPET002",30), rep("RS_fPET003",30),rep("RS_fPET004",30), rep("RS_fPET005",30),rep("RS_fPET006",30),rep("RS_fPET007",30),rep("RS_fPET008",30),rep("RS_fPET010",30),
rep("RS_fPET011",30),rep("RS_fPET012",30),rep("RS_fPET014",30),rep("RS_fPET015",30),rep("RS_fPET016",30),
rep("RS_fPET017",30),
rep("RS_fPET018",30),rep("RS_fPET019",30),rep("RS_fPET020",30), rep("RS_fPET021",30),rep("RS_fPET022",30),rep("RS_fPET023",30),rep("RS_fPET024",30),rep("RS_fPET025",30),
rep("RS_fPET026",30),rep("RS_fPET027",30),rep("RS_fPET028",30),rep("RS_fPET029",30))
#add record id and convert record id into factor
CogNETp$record_id<-subject_id
CogNETp$record_id<-factor(CogNETp$record_id)
#add diagnose to_id gm removed data set
CogNETp$diagnose[CogNETp$record_id=="RS_fPET001"|CogNETp$record_id=="RS_fPET002"|CogNETp$record_id=="RS_fPET003"|CogNETp$record_id=="RS_fPET005"|CogNETp$record_id=="RS_fPET007"|CogNETp$record_id=="RS_fPET008"|CogNETp$record_id=="RS_fPET010"|CogNETp$record_id=="RS_fPET014"|CogNETp$record_id=="RS_fPET016"|CogNETp$record_id=="RS_fPET022"|CogNETp$record_id=="RS_fPET025"|CogNETp$record_id=="RS_fPET026"|CogNETp$record_id=="RS_fPET028"|CogNETp$record_id=="RS_fPET029"]= "PD"
CogNETp$diagnose[CogNETp$record_id=="RS_fPET004"|CogNETp$record_id=="RS_fPET006"|CogNETp$record_id=="RS_fPET011"|CogNETp$record_id=="RS_fPET012"|CogNETp$record_id=="RS_fPET015"|CogNETp$record_id=="RS_fPET017"|CogNETp$record_id=="RS_fPET018"|CogNETp$record_id=="RS_fPET019"|CogNETp$record_id=="RS_fPET020"|CogNETp$record_id=="RS_fPET021"|CogNETp$record_id=="RS_fPET023"|CogNETp$record_id=="RS_fPET024"|CogNETp$record_id=="RS_fPET027"]= "HC"
#convert diagnose into factor
CogNETp$diagnose<-factor(CogNETp$diagnose)
#show record ids
CogNETp$record_id
## [1] RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001
## [7] RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001
## [13] RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001
## [19] RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001
## [25] RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001 RS_fPET001
## [31] RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002
## [37] RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002
## [43] RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002
## [49] RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002
## [55] RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002 RS_fPET002
## [61] RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003
## [67] RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003
## [73] RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003
## [79] RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003
## [85] RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003 RS_fPET003
## [91] RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004
## [97] RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004
## [103] RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004
## [109] RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004
## [115] RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004 RS_fPET004
## [121] RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005
## [127] RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005
## [133] RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005
## [139] RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005
## [145] RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005 RS_fPET005
## [151] RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006
## [157] RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006
## [163] RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006
## [169] RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006
## [175] RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006 RS_fPET006
## [181] RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007
## [187] RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007
## [193] RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007
## [199] RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007
## [205] RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007 RS_fPET007
## [211] RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008
## [217] RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008
## [223] RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008
## [229] RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008
## [235] RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008 RS_fPET008
## [241] RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010
## [247] RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010
## [253] RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010
## [259] RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010
## [265] RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010 RS_fPET010
## [271] RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011
## [277] RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011
## [283] RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011
## [289] RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011
## [295] RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011 RS_fPET011
## [301] RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012
## [307] RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012
## [313] RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012
## [319] RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012
## [325] RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012 RS_fPET012
## [331] RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014
## [337] RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014
## [343] RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014
## [349] RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014
## [355] RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014 RS_fPET014
## [361] RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015
## [367] RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015
## [373] RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015
## [379] RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015
## [385] RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015 RS_fPET015
## [391] RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016
## [397] RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016
## [403] RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016
## [409] RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016
## [415] RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016 RS_fPET016
## [421] RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017
## [427] RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017
## [433] RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017
## [439] RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017
## [445] RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017 RS_fPET017
## [451] RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018
## [457] RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018
## [463] RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018
## [469] RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018
## [475] RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018 RS_fPET018
## [481] RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019
## [487] RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019
## [493] RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019
## [499] RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019
## [505] RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019 RS_fPET019
## [511] RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020
## [517] RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020
## [523] RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020
## [529] RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020
## [535] RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020 RS_fPET020
## [541] RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021
## [547] RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021
## [553] RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021
## [559] RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021
## [565] RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021 RS_fPET021
## [571] RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022
## [577] RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022
## [583] RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022
## [589] RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022
## [595] RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022 RS_fPET022
## [601] RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023
## [607] RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023
## [613] RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023
## [619] RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023
## [625] RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023 RS_fPET023
## [631] RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024
## [637] RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024
## [643] RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024
## [649] RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024
## [655] RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024 RS_fPET024
## [661] RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025
## [667] RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025
## [673] RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025
## [679] RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025
## [685] RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025 RS_fPET025
## [691] RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026
## [697] RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026
## [703] RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026
## [709] RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026
## [715] RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026 RS_fPET026
## [721] RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027
## [727] RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027
## [733] RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027
## [739] RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027
## [745] RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027 RS_fPET027
## [751] RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028
## [757] RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028
## [763] RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028
## [769] RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028
## [775] RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028 RS_fPET028
## [781] RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029
## [787] RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029
## [793] RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029
## [799] RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029
## [805] RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029 RS_fPET029
## 27 Levels: RS_fPET001 RS_fPET002 RS_fPET003 RS_fPET004 ... RS_fPET029
#add timepoints
time<-c(rep(61:90, 27)) #repeat time points 61- 90 16 x
CogNETp$time<-time
#insert MCI based on ID according to RSfPET demographics script
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET003"]="PD-MCI"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET005"]="PD-MCI"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET010"]="PD-MCI"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET014"]="PD-MCI"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET025"]="PD-MCI"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET001"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET002"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET007"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET008"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET016"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET022"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET026"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET028"]="PD-NC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET029"]="PD-NC"
#only non-MCI controls
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET004"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET011"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET012"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET015"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET018"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET019"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET020"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET021"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET023"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET024"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET027"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET017"]="HC"
CogNETp$MCI_NC_HC[CogNETp$record_id =="RS_fPET006"]="HC"
#HC RSfPET017 and RSfPET06 have MCI
CogNETp$MCI_NC_HC<-factor(CogNETp$MCI_NC_HC)
for (i in 1:ncol(CogNETp[,1:20])) {
MeanTime <- aggregate(CogNETp[,i] ~ time + MCI_NC_HC, data = CogNETp, FUN = mean)
MeanTime_df <- cbind(MeanTime[-ncol(MeanTime)], MeanTime[[ncol(MeanTime)]])
colnames(MeanTime_df)<-c("time","diagnose",colnames(CogNETp)[i])
for (j in 3:ncol(MeanTime_df)){
plot <- ggplot(CogNETp, aes(x = time, y = CogNETp[,i])) +
geom_line(aes(color = MCI_NC_HC, group = record_id), linetype = 3,, linewidth = 1) +
geom_line(data = MeanTime_df, aes(x = time, y = MeanTime_df[,j], color=diagnose), linewidth = 2) +
scale_color_manual(values=c("#440154","#BFD200","#1F78B4"))+
theme_bw()+
labs(title = "", x = "Time series (Frames)",
y = "[18F]-FDG Uptake")+
theme(legend.title = element_blank(),legend.text = element_text(size=20),
plot.margin=margin(50,20,20,20))+
theme(axis.text=element_text(size=25),
axis.title=element_text(size=30,face="bold", margin=margin(t=20)))+
show(plot)
ggsave(plot,file=paste0("CogNETp",width=12, height=9,units="in", dpi=400,limitsize=FALSE,colnames(CogNETp)[i],".png"))}}
## function (x, y, ...)
## UseMethod("plot")
## <bytecode: 0x123535270>
## <environment: namespace:base>
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
#calculate variation coefficient
#calculate variation in subjects
Subjects<-c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029")
#define region names
region<-c("MPFC","LPL","LPR","PCC","SML","SMR","SMS","ACC","AIL","AIR","PFCL","PFCR","SMGL","SMGR","LPFCL","PPCL","LPFCR","PPCR","CA","CP")
#define matrix for variation coefficient
Vk<-matrix(ncol=27, nrow=20)
#loop over subjects to fill matrix with variataion coefficient per subject
for(i in 1:length(Subjects)){
Vk[1,i]<- sd(subset(CogNETp,record_id==Subjects[i])$MPFC)/mean(subset(CogNETp,record_id==Subjects[i])$MPFC)*100
Vk[2,i]<- sd(subset(CogNETp,record_id==Subjects[i])$LPL)/mean(subset(CogNETp,record_id==Subjects[i])$LPL)*100
Vk[3,i]<- sd(subset(CogNETp,record_id==Subjects[i])$LPR)/mean(subset(CogNETp,record_id==Subjects[i])$LPR)*100
Vk[4,i]<- sd(subset(CogNETp,record_id==Subjects[i])$PCC)/mean(subset(CogNETp,record_id==Subjects[i])$PCC)*100
Vk[5,i]<- sd(subset(CogNETp,record_id==Subjects[i])$SML)/mean(subset(CogNETp,record_id==Subjects[i])$SML)*100
Vk[6,i]<- sd(subset(CogNETp,record_id==Subjects[i])$SMR)/mean(subset(CogNETp,record_id==Subjects[i])$SMR)*100
Vk[7,i]<- sd(subset(CogNETp,record_id==Subjects[i])$SMS)/mean(subset(CogNETp,record_id==Subjects[i])$SMS)*100
Vk[8,i]<- sd(subset(CogNETp,record_id==Subjects[i])$ACC)/mean(subset(CogNETp,record_id==Subjects[i])$ACC)*100
Vk[9,i]<- sd(subset(CogNETp,record_id==Subjects[i])$AIL)/mean(subset(CogNETp,record_id==Subjects[i])$AIL)*100
Vk[10,i]<- sd(subset(CogNETp,record_id==Subjects[i])$AIR)/mean(subset(CogNETp,record_id==Subjects[i])$AIR)*100
Vk[11,i]<- sd(subset(CogNETp,record_id==Subjects[i])$PFCL)/mean(subset(CogNETp,record_id==Subjects[i])$PFCL)*100
Vk[12,i]<- sd(subset(CogNETp,record_id==Subjects[i])$PFCR)/mean(subset(CogNETp,record_id==Subjects[i])$PFCR)*100
Vk[13,i]<- sd(subset(CogNETp,record_id==Subjects[i])$SMGL)/mean(subset(CogNETp,record_id==Subjects[i])$SMGL)*100
Vk[14,i]<- sd(subset(CogNETp,record_id==Subjects[i])$SMGR)/mean(subset(CogNETp,record_id==Subjects[i])$SMGR)*100
Vk[15,i]<- sd(subset(CogNETp,record_id==Subjects[i])$LPFCL)/mean(subset(CogNETp,record_id==Subjects[i])$LPFCL)*100
Vk[16,i]<- sd(subset(CogNETp,record_id==Subjects[i])$PPCL)/mean(subset(CogNETp,record_id==Subjects[i])$PPCL)*100
Vk[17,i]<- sd(subset(CogNETp,record_id==Subjects[i])$LPFCR)/mean(subset(CogNETp,record_id==Subjects[i])$LPFCR)*100
Vk[18,i]<- sd(subset(CogNETp,record_id==Subjects[i])$PPCR)/mean(subset(CogNETp,record_id==Subjects[i])$PPCR)*100
Vk[19,i]<- sd(subset(CogNETp,record_id==Subjects[i])$CA)/mean(subset(CogNETp,record_id==Subjects[i])$CA)*100
Vk[20,i]<- sd(subset(CogNETp,record_id==Subjects[i])$CP)/mean(subset(CogNETp,record_id==Subjects[i])$CP)*100
}
#transform into dataframe
Vk_data<-data.frame(
Vk)
colnames(Vk_data)<-Subjects
tVk_data<-t(Vk_data)
colnames(tVk_data)<-region
#in groups HC vs PD
ID_HCPD<-c("PD","PD","PD","HC","PD","HC","PD","PD","PD","HC","HC","PD","HC","PD","HC","HC","HC","HC","HC","PD","HC","HC","PD","PD","HC","PD","PD")
tVk_group<-as.data.frame(tVk_data)
tVk_group<-data.frame(lapply(tVk_group,as.numeric))
tVk_group$ID_HCPD<-factor(ID_HCPD)
summary(tVk_group$ID_HCPD)
## HC PD
## 13 14
variables<- colnames(tVk_group[,1:12])
grp<-levels(tVk_group$ID_HCPD)
#function
perm_test<- function(variable, group_var){
group_var<-as.factor(group_var)
grp_levels<-levels(group_var)
obs_diff<-tapply(variable,group_var,mean, na.rm=TRUE)
obs_diff_A<-obs_diff[grp_levels[2]]-obs_diff[grp_levels[1]]
perm_diffsA<-replicate(10000,{
permuted_groups<-sample(group_var)
perm_mean_diff_A<-tapply(variable,permuted_groups, mean, na.rm=TRUE)
perm_mean_diff_A[grp_levels[2]]-perm_mean_diff_A[grp_levels[1]]
})
print(length(perm_diffsA))
#p-values
p_value_A<-mean(abs(perm_diffsA)>=abs(obs_diff_A))
return(c(p_value_A))
}
#save results
results<-data.frame(Variable=character(), p_value=numeric(), stringsAsFactors = FALSE)
#loop over variables
for (var in variables){
p_values<-perm_test(tVk_group[[var]], tVk_group$ID_HCPD)
results<-rbind(results, data.frame(Variable=var, p_value_A=p_values[1]))
}
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
print(results)
## Variable p_value_A
## 1 MPFC 0.8334
## 2 LPL 0.4866
## 3 LPR 0.2651
## 4 PCC 0.0550
## 5 SML 0.9459
## 6 SMR 0.2244
## 7 SMS 0.0358
## 8 ACC 0.3419
## 9 AIL 0.9251
## 10 AIR 0.6611
## 11 PFCL 0.8232
## 12 PFCR 0.8606
#for multple boxplots we need additional long format
ggplot(long, aes(x=region, y=value, fill=ID_HCPD)) +
geom_half_violin(side="r",trim=FALSE, alpha=0.5, show.legend=FALSE)+
geom_half_boxplot(side="l")+
facet_wrap(~region, scale="free")+
# coord_cartesian(ylim=c(min(long$value), max(long$value)+1))+
scale_fill_manual(labels=c("HC","PD"),values=c("#1F78B4", "#BFD200"))+
theme_minimal()+
labs(title = "", x = "Region",
y = "VC")+
theme(text = element_text(size = 15),
title = element_text(size = 20),
axis.text.x=element_blank(),
axis.text.y = element_text(size=9),
legend.title=element_blank())
## Error in eval(expr, envir, enclos): object 'long' not found
ggsave("Vk_atlas_boxplots.png")
## Saving 7 x 5 in image
## Error in `geom_line()`:
## ! Problem while computing aesthetics.
## ℹ Error occurred in the 1st layer.
## Caused by error in `[.data.frame`:
## ! undefined columns selected
#with Cog Z from demographic script
load("RSfPET.RData")
#in whole gruop with cog z
plot(tVk_group$SMS,subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$Cog_z)
model<-lm(subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$Cog_z~tVk_group$SMS)
abline(model)
summary(model)
##
## Call:
## lm(formula = subset(RSfPET, record_id %in% c("RS_fPET001", "RS_fPET002",
## "RS_fPET003", "RS_fPET004", "RS_fPET005", "RS_fPET006", "RS_fPET007",
## "RS_fPET008", "RS_fPET010", "RS_fPET011", "RS_fPET012", "RS_fPET014",
## "RS_fPET015", "RS_fPET016", "RS_fPET017", "RS_fPET018", "RS_fPET019",
## "RS_fPET020", "RS_fPET021", "RS_fPET022", "RS_fPET023", "RS_fPET024",
## "RS_fPET025", "RS_fPET026", "RS_fPET027", "RS_fPET028", "RS_fPET029"))$Cog_z ~
## tVk_group$SMS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.77102 -0.35105 0.00295 0.48536 1.92220
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1496 0.6271 1.833 0.0787 .
## tVk_group$SMS -0.7714 0.3381 -2.281 0.0313 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8706 on 25 degrees of freedom
## Multiple R-squared: 0.1723, Adjusted R-squared: 0.1392
## F-statistic: 5.205 on 1 and 25 DF, p-value: 0.03131
#in whole group with moca
plot(tVk_group$SMS,subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$MoCa_Sum)
model<-lm(subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$MoCa_Sum~tVk_group$SMS)
abline(model)
summary(model)
##
## Call:
## lm(formula = subset(RSfPET, record_id %in% c("RS_fPET001", "RS_fPET002",
## "RS_fPET003", "RS_fPET004", "RS_fPET005", "RS_fPET006", "RS_fPET007",
## "RS_fPET008", "RS_fPET010", "RS_fPET011", "RS_fPET012", "RS_fPET014",
## "RS_fPET015", "RS_fPET016", "RS_fPET017", "RS_fPET018", "RS_fPET019",
## "RS_fPET020", "RS_fPET021", "RS_fPET022", "RS_fPET023", "RS_fPET024",
## "RS_fPET025", "RS_fPET026", "RS_fPET027", "RS_fPET028", "RS_fPET029"))$MoCa_Sum ~
## tVk_group$SMS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.1125 -1.1603 -0.0165 2.0909 4.8614
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.912 2.442 13.475 5.76e-13 ***
## tVk_group$SMS -3.370 1.317 -2.559 0.0169 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.391 on 25 degrees of freedom
## Multiple R-squared: 0.2075, Adjusted R-squared: 0.1759
## F-statistic: 6.548 on 1 and 25 DF, p-value: 0.01694
#whole group Cogz scores (demographics file)
ggplot(data=tVk_group,aes(x=SMS, y=subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$Cog_z, color=ID_HCPD, shape=RSfPET$MCI_NC_HC)) +
geom_point(size=3)+
geom_smooth(data=tVk_group,aes(group=1),method='lm', formula= y~x,show.legend = FALSE, color="black", se=FALSE)+
scale_color_manual(labels=c("HC","PD"),values=c("#1F78B4", "#BFD200"))+
theme_minimal()+
theme(legend.position = "none", legend.title=element_blank())+
labs(title = "", x = "SMS VC",
y = "Cognition z-score")+
theme(text = element_text(size = 30),
axis.text.x = element_text(size=20),
axis.text.y = element_text(size=20))+
scale_shape_manual(values=c("HC"=16, "MCI"=17, "PD-NC"=16,"PD-MCI"=17))+
geom_text(x=1.5,y=-1.0, label="r = -0.42, p = 0.031", color="black", size=7)
## Warning: The following aesthetics were dropped during statistical transformation: shape.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("Vk_SMS_CogZ.png")
## Saving 7 x 5 in image
## Warning: The following aesthetics were dropped during statistical transformation: shape.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
cor.test(subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$Cog_z,tVk_group$SMS, method="pearson")
##
## Pearson's product-moment correlation
##
## data: subset(RSfPET, record_id %in% c("RS_fPET001", "RS_fPET002", "RS_fPET003", "RS_fPET004", "RS_fPET005", "RS_fPET006", "RS_fPET007", "RS_fPET008", "RS_fPET010", "RS_fPET011", "RS_fPET012", "RS_fPET014", "RS_fPET015", "RS_fPET016", "RS_fPET017", "RS_fPET018", "RS_fPET019", "RS_fPET020", "RS_fPET021", "RS_fPET022", "RS_fPET023", "RS_fPET024", "RS_fPET025", "RS_fPET026", "RS_fPET027", "RS_fPET028", "RS_fPET029"))$Cog_z and tVk_group$SMS
## t = -2.2814, df = 25, p-value = 0.03131
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.68678067 -0.04166061
## sample estimates:
## cor
## -0.4151028
cor.test(subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$MoCa_Sum,tVk_group$SMS, method="spearman")
## Warning in cor.test.default(subset(RSfPET, record_id %in% c("RS_fPET001", :
## Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: subset(RSfPET, record_id %in% c("RS_fPET001", "RS_fPET002", "RS_fPET003", "RS_fPET004", "RS_fPET005", "RS_fPET006", "RS_fPET007", "RS_fPET008", "RS_fPET010", "RS_fPET011", "RS_fPET012", "RS_fPET014", "RS_fPET015", "RS_fPET016", "RS_fPET017", "RS_fPET018", "RS_fPET019", "RS_fPET020", "RS_fPET021", "RS_fPET022", "RS_fPET023", "RS_fPET024", "RS_fPET025", "RS_fPET026", "RS_fPET027", "RS_fPET028", "RS_fPET029"))$MoCa_Sum and tVk_group$SMS
## S = 4576.6, p-value = 0.04033
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.397